**Supplementary Analysis & Results**
version
# Dace ApŔvalka @MRC-CBU, 2019-2020
# ----------------------------------------------------------------------
# LIBRARIES
# ----------------------------------------------------------------------
shhh <- suppressPackageStartupMessages # It's a library, so shhh!
shhh(suppressWarnings(library(Hmisc))) # for several stats things
shhh(suppressWarnings(library(dplyr))) # for data manipulation
shhh(suppressWarnings(library(plyr))) # for ddply (aggregating SST results)
shhh(suppressWarnings(library(pastecs))) # for summary statistic
shhh(suppressWarnings(library(ggplot2))) # for plotting results
shhh(suppressWarnings(library(cowplot))) # for adding plots together and setting different widths
shhh(suppressWarnings(library(repr))) # for changing plot size (default too large in notebook)
shhh(suppressWarnings(library(gridExtra))) # for adding plots together
shhh(suppressWarnings(library(grid))) # for plot grid title formatting
shhh(suppressWarnings(library(aplpack))) # for bi-variate outlier, bgplot
shhh(suppressWarnings(library(robust))) # for robust correlation
shhh(suppressWarnings(library(tidyr))) # for creating long table (function gather)
shhh(suppressWarnings(library(gtools))) # for converting pvalues to stars
shhh(suppressWarnings(library(ggpubr))) # for adding statistics to plots
shhh(suppressWarnings(library(rstatix))) # for Basic Statistical Tests
shhh(suppressWarnings(library(Rmisc))) # for getting summary data frame
shhh(suppressWarnings(library(WRS2))) # for percentage bend correlation coefficient
# https://link-springer-com.ezp.lib.cam.ac.uk/content/pdf/10.1007/BF02294395.pdf
shhh(suppressWarnings(library(WRS))) # for Wilcox correlation functions (skipped correlations)
# Wilcox R. R. (2005). Introduction to Robust Estimation and Hypothesis Testing.
# New York, NY: Elsevier Academic Press
# https://www-ncbi-nlm-nih-gov.ezp.lib.cam.ac.uk/pmc/articles/PMC3342588/
shhh(suppressWarnings(library(MASS))) # WRS2 uses it. Its select function conflicts with dplyr select!
options(warn=-1)
# STOP-SIGNAL ANALYSIS
# Based on: Verbruggen et al (2019). https://doi.org/10.7554/eLife.46323.001
# ------------------------------------
funcSignal <- function(data) {
# all trials with a Response (any)
data$resp <- ifelse(data$gRT == "0", 0, 1)
# Stop-signal trials
# ----------------------------------------------------------------------
# All
signal <- subset(data, gkind == '1')
# probability of responding to the Stop-Signal
presp <- mean(signal$resp)
# mean Stop-Signal Delay (of all Stop-Signal trials)
ssd <- mean(signal$gintVal)
# unsuccessful Stop RTs
signal.resp <- subset(signal, resp == '1')
signal.resp.rt <- mean(signal.resp$gRT)
# Go trials
# ----------------------------------------------------------------------
# All
nosignal <- subset(data, gkind == '0')
# With response
nosignal_resp <- subset(nosignal, resp == '1')
# probability of Go omossions
pmiss <- 1-mean(nosignal$resp)
# Correct Go trials
nosignal_resp$error <-
ifelse(nosignal_resp$gtheButton == nosignal_resp$gresp, 0, 1)
# Probability of Go choice errors
perror <- mean(nosignal_resp$error)
# for the missed responses, set RT to max RT of the subject
nosignal$gRT <-
ifelse(nosignal$gRT == 0, max(nosignal_resp$gRT), nosignal$gRT)
# RT for Go trials with a response
nosignal.resp.rt <- mean(nosignal_resp$gRT)
# RT for correct Go trials
nosignal.correctresp <- subset(nosignal_resp, error == '0')
nosignal.correctresp.rt <- mean(nosignal.correctresp$gRT)
# SSRT
# ----------------------------------------------------------------------
# All Go trials are included when the nth RT is determined
## determine nth RT
nthRT_all <- quantile(nosignal$gRT, probs = presp, type = 6)
## SSRT(integration) = nthRT - ssd
SSRTall <- nthRT_all - ssd
# Race Check
# ----------------------------------------------------------------------
# All Go trials
# "..this comparison should include all trials with a response
# (including choice errors and premature responses).."
# (Verbruggen et al., 2019))
race.check <- nosignal.resp.rt - signal.resp.rt
# Output
# ----------------------------------------------------------------------
return(
data.frame(
Ntrials = nrow(data),
NStop = sum(signal$gkind),
ssPerc = sum(signal$gkind)/nrow(data)*100,
presp, pmiss, perror,
SSD = ssd,
SSRT = SSRTall,
sRT = signal.resp.rt,
goRTcorr = nosignal.correctresp.rt,
goRTall = nosignal.resp.rt,
raceCheck = race.check
)
)
}
# PLOT MEANS, Group mean (with 95%CI) and subject values
# =======================================================
plotResults <- function(dataset, x, group, title, units, extraTxt) {
ggplot(dataset, aes(group, x)) +
geom_point(
colour = "black", alpha = .9, fill = 'grey',
size = 1, stroke = 0.1, shape = 21,
position = position_jitter(height = 0, width = 0.1)
) +
stat_summary(
fun.data = mean_cl_boot,
geom = "pointrange",
fill = 'red', alpha = 0.8, color = "black",
size = 0.6, stroke = 0.3, shape = 21,
position = position_nudge(x = -0.2)
) +
# scale_x_discrete() +
labs(x = extraTxt, y = units) + ggtitle(title) +
theme_minimal() +
theme(text = element_text(size = 9),
plot.title = element_text(hjust = 0.5, size = 9, face="bold"),
plot.background = element_rect(color = "lightgrey", size = 0.4, linetype = "dotted"),
panel.grid.major = element_blank(), panel.grid.minor = element_blank()
)
}
# GET OUTLIERS (from a variable pair, for correlations)
# =======================================================
get_outliers <-
function(var1,var1name,var1y,
var2,var2name,var2y,
disp) {
par(mfrow = c(1, 3)) # to put the 2 univariate and 1 bivariate plots together
# UNIVARIATE, boxplot method
# ------------------------------------
bpVar1 <- boxplot(var1,main = var1name,ylab = var1y,plot = disp)
bpVar2 <- boxplot(var2,main = var2name,ylab = var2y,plot = disp)
# oVar1
ifelse(length(bpVar1$out) == 0,
oVar1 <- vector(mode = "numeric", length = 0),
oVar1 <- which(var1 %in% bpVar1$out)
)
# oVar2
ifelse(length(bpVar2$out) == 0,
oVar2 <- vector(mode = "numeric", length = 0),
oVar2 <- which(var2 %in% bpVar2$out)
)
# BI-VARIATE, bagplot method
# ------------------------------------
# A bagplot is a bivariate generalization of the boxplot. A bivariate boxplot.
bpl <- compute.bagplot(var1, var2)
# oBiv
ifelse(is.null(bpl$pxy.outlier),
oBiv <- vector(mode = "numeric", length = 0),
oBiv <- which(var1 %in% bpl$pxy.outlier[, 1])
)
bp <-
bagplot(
var1,var2,
ylab = var2name,xlab = var1name,
main = 'bi-variate',
cex = 2,
create.plot = disp
)
# Report results
# ------------------------------------
if(disp){
# var1
print(sprintf('%s outliers: %d', var1name, length(bpVar1$out)))
print(oVar1)
# var2
print(sprintf('%s outliers: %d', var2name, length(bpVar2$out)))
print(oVar2)
# bi-var
print(sprintf('Bi-variate outliers: %d', length(oBiv)))
print(oBiv)
}
# Return results
# ------------------------------------
all <- c(oVar1, oVar2, oBiv)
univariate <- c(oVar1, oVar2)
bivariate <- oBiv
return(list(all, univariate, bivariate))
}
# DO CORRELATIONS
# Depending on the data, 3 types of correlations possible. Following recommendations by Pernet et al.(2013)
# =======================================================
doCorrelation <- function(var1, var2, outliers) {
# WHICH CORRELATION
# ------------------------------------
isOutliers <- length(outliers[[1]]) > 0
isUnivariate <- length(outliers[[2]]) > 0
isBivariate <- length(outliers[[3]]) > 0
isNormal <- (shapiro.test(var1)$p.value > 0.05 & shapiro.test(var2)$p.value > 0.05)
# if is Bivariate do Spearman skipped, using the minimum covariance determinant (MCD) estimator
if (isBivariate) {
# corRes <-
# scorci(var1,var2,nboot = 1000,corfun = spear,plotit = FALSE)
# r <- corRes$cor.est
# p <- corRes$p.value
# f <- 'Spearman skipped'
corRes <-
mscor(data.frame(var1,var2), corfun = spear, cop = 3)
r <- corRes$cor[2]
p <- 2*pt(-abs(corRes$test.stat[2]), df = length(var1)-1)
f <- 'Spearman skipped'
subs <- 'ss'
}
# if is not Bivar but is Univar or is not Normal, do 20% Bend
if (!isBivariate & (isUnivariate | !isNormal)) {
corRes <-
pbcor(var1,var2,beta = 0.2)
r <- corRes$cor
p <- corRes$p.value
f <- 'Percentage-bend'
subs <- 'pb'
}
# if no outliers and is normal do Pearson
if (!isOutliers & isNormal) {
corRes <-
rcorr(var1,var2)
r <- corRes$r[2]
p <- corRes$P[2]
f <- 'Pearson'
subs <- ''
}
pval <- ifelse(p < 0.001, "p < 0.001", sprintf("p = %.3f", p))
corResTxt <- bquote(.(f) ~ "correlation" ~ r[.(subs)] == .(sprintf("%.3f, %s", r, pval)))
return(corResTxt)
}
# PLOT CORRELATIONS (with 95%CI)
# ------------------------------------
plotCorrelation <- function(var1, var2, var1name, var2name, resTXT,
pointsize, txtsize, outliers = NULL, plotoutliers = FALSE) {
if(!is.null(outliers) & length(outliers[[1]]) > 0) {
out1 <- var1[outliers[[1]]]
out2 <- var2[outliers[[1]]]
outdata <- data.frame(out1, out2)
var1 <- var1[-c(outliers[[1]])]
var2 <- var2[-c(outliers[[1]])]
ifelse(plotoutliers == TRUE,
addTxt <- sprintf('\n Outliers (n=%d) displayed in red',length(unique(outliers[[1]]))),
addTxt <- sprintf('\n Outliers (n=%d) not displayed',length(unique(outliers[[1]])))
)
resTXT <- bquote(atop(.(resTXT), .(addTxt)))
#resTXT <- paste(resTXT, sprintf('\n Outliers (n=%d) not displayed',length(unique(outliers[[1]]))))
}
dataset <- data.frame(var1, var2)
corplot <- ggplot(dataset, aes(var1, var2)) +
geom_smooth(
method = lm, formula = y ~ x, level = 0.95,
color = "black", fill = "grey", size = 0.2
) +
geom_point(
colour = "black", alpha = .9, fill = 'grey',
size = pointsize, stroke = 0.2, shape = 21
) +
labs(x = var1name, y = var2name, title = resTXT) +
theme_minimal() +
theme(text = element_text(size = txtsize),
plot.title = element_text(hjust = 0.5, size = 11, face="bold"),
plot.background = element_rect(colour = "lightgrey", size = 0.4, linetype = "dotted"),
panel.grid.major = element_blank(), panel.grid.minor = element_blank()
)
if(!is.null(outliers) & length(outliers[[1]]) > 0 & plotoutliers == TRUE) {
corplot <- corplot + geom_point(data = outdata, aes(out1, out2), color = 'red')
}
return(corplot)
}
# 2 by 2 Repeated measures ANOVA
# =======================================================
# using my function described here https://dcdace.net/2020/10/12/rm_anova/
source("https://raw.githubusercontent.com/dcdace/R_functions/main/rm_2by2_anova.R")
# Suppress function output
# =======================================================
# I got it from here https://stackoverflow.com/questions/2723034/suppress-output-of-a-function
hush <- function(code){
sink("NUL") # use /dev/null in UNIX
tmp = code
sink()
return(tmp)
}
# ----------------------------------------------------------------------
# LOAD DATA
# ----------------------------------------------------------------------
# SST
dataSST <- read.csv('https://raw.githubusercontent.com/dcdace/tmp/master/data/SSTraw.csv')
# exclude Go only runs from the analysis
dataSST <- subset(dataSST, cond == '2')
# TNT
dataTNT <- read.csv('https://raw.githubusercontent.com/dcdace/tmp/master/data/TNT.csv')
# turn test values into percentages
dataTNT[, 10:ncol(dataTNT)] <- dataTNT[, 10:ncol(dataTNT)] * 100
# tfSSRT: SSRTs with Stop trigger failures
data_tfSSRT <- read.csv('https://raw.githubusercontent.com/dcdace/tmp/master/data/TF_SSRT.csv')
# ======================================================================
# SST RESULTS
# ======================================================================
# By Subject aggregated
# ----------------------------------------------------------------------
dataSSTsubj <- ddply(dataSST, .(sNR), funcSignal)
# display
options(width = 118)
print(round(dataSSTsubj,2))
# Descriptives
# ----------------------------------------------------------------------
dscrSST <- stat.desc(dataSSTsubj[, 2:ncol(dataSSTsubj)], basic = TRUE, desc = TRUE, p = 0.95)
print(round(dscrSST, 2))
# Plot results
# ----------------------------------------------------------------------
options(repr.plot.width = 10, repr.plot.height = 10, repr.plot.res = 200)
# 1. Probability of Go omissions (no response)
p1 <- plotResults(dataSSTsubj, dataSSTsubj$pmiss, 1, "p(miss|Go)", "p",
sprintf("M (SD) = %.3f (%.2f), Mdn = %.2f",
dscrSST["mean","pmiss"], dscrSST["std.dev","pmiss"], dscrSST["median","pmiss"])) +
scale_x_discrete()
# 2. Probability of choice errors on Go trials
p2 <- plotResults(dataSSTsubj, dataSSTsubj$perror, 1, "p(error|Go)", "p",
sprintf("M (SD) = %.2f (%.2f), Mdn = %.2f",
dscrSST["mean","perror"], dscrSST["std.dev","perror"],
dscrSST["median","perror"])) +
scale_x_discrete()
# 3. Probability of responding on a Stop trial
# recommended to refrain from estimating individual SSRTs when p (respond|signal)
# is lower than 0.25 or higher than 0.75
p3 <- plotResults(dataSSTsubj, dataSSTsubj$presp, 1, "p(respond|signal)", "p",
sprintf("M (SD) = %.2f (%.2f), Mdn = %.2f\nmin = %.2f, max = %.2f",
dscrSST["mean","presp"], dscrSST["std.dev","presp"], dscrSST["median","presp"],
dscrSST["min","presp"], dscrSST["max","presp"])) +
scale_x_discrete()
# 4. Average Stop-Signal Delay
p4 <- plotResults(dataSSTsubj, dataSSTsubj$SSD, 1, "SSD", "ms",
sprintf("M (SD) = %.2f (%.2f) ms, Mdn = %.2f\nmin = %.2f, max = %.2f",
dscrSST["mean","SSD"], dscrSST["std.dev","SSD"], dscrSST["median","SSD"],
dscrSST["min","SSD"], dscrSST["max","SSD"])) +
scale_x_discrete()
# 5. Mean RT on correct Go trials
p5 <- plotResults(dataSSTsubj, dataSSTsubj$goRTcorr, 1, "correct Go RT", "ms",
sprintf("M (SD) = %.2f (%.2f) ms",
dscrSST["mean","goRTcorr"], dscrSST["std.dev","goRTcorr"])) +
scale_x_discrete()
# 6. RT of unsuccessful Stop trials
p6 <- plotResults(dataSSTsubj, dataSSTsubj$sRT, 1, "failed Stop RT", "ms",
sprintf("M (SD) = %.2f (%.2f) ms",
dscrSST["mean","sRT"], dscrSST["std.dev","sRT"])) +
scale_x_discrete()
# 7. Race Check (with correct Go)
p7 <- plotResults(dataSSTsubj, dataSSTsubj$raceCheck, 1, "Race check", "ms",
sprintf("M (SD) = %.2f (%.2f) ms",
dscrSST["mean","raceCheck"], dscrSST["std.dev","raceCheck"])) +
scale_x_discrete() + geom_hline(aes(yintercept = 0), size = 0.1)
# 8. Stop-signal reaction time (SSRT)
# integration method (with replacement of Go omissions) when
p8 <- plotResults(dataSSTsubj, dataSSTsubj$SSRT, 1, "SSRT", "ms",
sprintf("M (SD) = %.2f (%.2f) ms",
dscrSST["mean","SSRT"], dscrSST["std.dev","SSRT"])) +
scale_x_discrete()
# 9. Change in Go RT across runs
dataSSTsubjRun <- ddply(dataSST, .(sNR, run), funcSignal)
p9.res.lm <- summary(lm(goRTall ~ run, data = dataSSTsubjRun))
p9.res.aov <- unlist(summary(aov(goRTall ~ run, data = dataSSTsubjRun)))
p9.rxt <- sprintf('B = %.3f, F = %.2f, p = %.3f',
p9.res.lm$coefficients["run", "Estimate"],
p9.res.aov["F value1"],
p9.res.lm$coefficients["run", "Pr(>|t|)"])
p9 <- ggplot(dataSSTsubjRun, aes(x = run, y = goRTall)) +
# subject lines
geom_smooth(method = lm, formula = y ~ x, se = FALSE, color = "black", size = 0.05) +
aes(fill = factor(sNR)) +
# group line
geom_smooth(colour = "black", method = lm, formula = y ~ x, size = 0.5, alpha = .5, fill = "black") +
# run means
stat_summary(fun.data = mean_cl_boot, geom = "pointrange", fill = 'red', alpha = 0.8,
color = "black", size = 0.6, stroke = 0.3, shape = 21) +
labs(x = "Run", y = "Go RT (ms)") + ggtitle("Go RT (all) change across runs") +
scale_x_continuous(breaks = unique(dataSSTsubjRun$run)) +
annotate("text", label= p9.rxt, x = 4, y = max(dataSSTsubjRun$goRTall), size = 3) +
theme_minimal() +
theme(text = element_text(size = 9),
plot.title = element_text(hjust = 0.5, size = 9, face="bold"),
plot.background = element_rect(color = "lightgrey", size = 0.4, linetype = "dotted"),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="none")
# 10. Inter-subject variability of SSD
dataSignal <- subset(dataSST, gkind == '1')
p10 <- plotResults(dataSignal, dataSignal$gintVal, dataSignal$sNR, "Inter-subject SSD", "ms", "") +
geom_hline(aes(yintercept = 0), size = 0.1) + scale_x_continuous(breaks = unique(dataSignal$sNR))
# 11. Inter-subject variability of RT on correct Go trials
#options(repr.plot.width = 8, repr.plot.height = 3)
dataNoSignal <- subset(dataSST, gkind == '0' & gRT != '0')
dataNoSignal$error <- ifelse(dataNoSignal$gtheButton == dataNoSignal$gresp, 0, 1)
dataNoSignalCorr <- subset(dataNoSignal, error == '0')
dataNoSignalCorr.rt <- mean(dataNoSignalCorr$gRT)
p11 <- plotResults(dataNoSignalCorr, dataNoSignalCorr$gRT, dataNoSignalCorr$sNR,
"Inter-subject correct Go RT", "ms", "") +
scale_x_continuous(breaks = unique(dataSignal$sNR))
## Put all plots in a grid and display
ly1 <- rbind(c(1,2,3),c(4,5,6), c(7,8,9), c(10,10,10), c(11,11,11))
grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, layout_matrix = ly1)
# ==================================================================
# TNT RESULTS
# ==================================================================
# display
# sp - Same probe test
# ip - Independent probe test
# uc - unconditionalised
# c - conditionalised
options(width = 114)
print(cbind(dataTNT[, 4:10], round(dataTNT[, 11:ncol(dataTNT)], 0)))
# Descriptives
# ----------------------------------------------------------------------
dscrTNT <- stat.desc(dataTNT[, 7:ncol(dataTNT)], basic = TRUE, desc = TRUE, p = 0.95)
options(width = 110)
print(round(dscrTNT, 0))
# Plot results
# ----------------------------------------------------------------------
options(repr.plot.width = 9, repr.plot.height = 6, repr.plot.res = 200)
# 1. Unconditionalised SP SIF
p12 <- plotResults(dataTNT, dataTNT$spSIFuc, 1, "SIF (SP)", "%",
sprintf("M (SD) = %.0f (%.0f) %%",
dscrTNT["mean","spSIFuc"], dscrTNT["std.dev","spSIFuc"])) +
scale_x_discrete() + geom_hline(aes(yintercept = 0), size = 0.1)
# 2. Unconditionalised IP SIF
p13 <- plotResults(dataTNT, dataTNT$ipSIFuc, 1, "SIF (IP)", "%",
sprintf("M (SD) = %.0f (%.0f) %%",
dscrTNT["mean","ipSIFuc"], dscrTNT["std.dev","ipSIFuc"])) +
scale_x_discrete() + geom_hline(aes(yintercept = 0), size = 0.1)
# 3. Unconditionalised SPIP SIF SP and IP
p14 <- plotResults(dataTNT, dataTNT$spipSIFuc, 1, "SIF (SP and IP)", "%",
sprintf("M (SD) = %.0f (%.0f) %%",
dscrTNT["mean","spipSIFuc"], dscrTNT["std.dev","spipSIFuc"])) +
scale_x_discrete() + geom_hline(aes(yintercept = 0), size = 0.1)
# 4. Conditionalised SP SIF
p15 <- plotResults(dataTNT, dataTNT$spSIFc, 1, "SIF (SP)", "%",
sprintf("M (SD) = %.0f (%.0f) %%",
dscrTNT["mean","spSIFc"], dscrTNT["std.dev","spSIFc"])) +
scale_x_discrete() + geom_hline(aes(yintercept = 0), size = 0.1)
# 5. Conditionalised IP SIF
p16 <- plotResults(dataTNT, dataTNT$ipSIFc, 1, "SIF (IP)", "%",
sprintf("M (SD) = %.0f (%.0f) %%",
dscrTNT["mean","ipSIFc"], dscrTNT["std.dev","ipSIFc"])) +
scale_x_discrete() + geom_hline(aes(yintercept = 0), size = 0.1)
# 6. Conditionalised SIF SP and IP
p17 <- plotResults(dataTNT, dataTNT$spipSIFc, 1, "SIF (SP and IP)", "%",
sprintf("M (SD) = %.0f (%.0f) %%",
dscrTNT["mean","spipSIFc"], dscrTNT["std.dev","spipSIFc"])) +
scale_x_discrete() + geom_hline(aes(yintercept = 0), size = 0.1) + theme(
panel.background = element_rect(fill = "lightyellow"))
# 7. Relationship between conditionalised and unconditionalised SIFs
outliers <- get_outliers(dataTNT$spipSIFuc, "unc. SIF", "%", dataTNT$spipSIFc, "cond. SIF", "", disp = FALSE)
corRes <- doCorrelation(dataTNT$spipSIFuc, dataTNT$spipSIFc, outliers)
p18 <- plotCorrelation(dataTNT$spipSIFuc,
dataTNT$spipSIFc,"Uncond. SIF","Cond. SIF",corRes, 1.8, 11, outliers) +
geom_vline(aes(xintercept = 0), size = 0.1) + geom_hline(aes(yintercept = 0), size = 0.1)
# 8. % of learned (included in analysis) items
p19 <- plotResults(dataTNT, dataTNT$criterion, 1, "% of learned items", "%",
sprintf("M (SD) = %.0f (%.0f) %%\nmin = %.0f, max = %.0f",
dscrTNT["mean","criterion"], dscrTNT["std.dev","criterion"],
dscrTNT["min","criterion"], dscrTNT["max","criterion"])) +
scale_x_discrete()
## Put all plots in a grid and display
grid.arrange(arrangeGrob(p12, p13, p14,
top = textGrob("Unconditionalised", gp=gpar(fontface="bold")), nrow = 1),
arrangeGrob(p15, p16, p17,
top = textGrob("Conditionalised", gp=gpar(fontface="bold")), nrow = 1),
arrangeGrob(p18, p19,
top="", layout_matrix = rbind(c(1,1,1,NA,2,2))))
# SIF effect (conditionalised, SP and IP)
# ----------------------------------------------------------------------
# http://www.sthda.com/english/wiki/one-sample-t-test-in-r
# Since the sample size is not large enough (less than 30, central limit theorem), we need to check whether
# the data follow a normal distribution.
# It's possible to use the Shapiro-Wilk normality test and to look at the normality plot.
#
# Shapiro-Wilk test:
# Null hypothesis: the data are normally distributed
# Alternative hypothesis: the data are not normally distributed
shapiro.test(dataTNT$spipSIFc)
# SIF data are normally distributed! Can do t.test. Otherwise would need to to wilcox.test
t.SIFc <- t.test(dataTNT$spipSIFc, mu = 0, alternative = "greater")
t.SIFc
sprintf("SIF effect: t(%d) = %.2f, p = %.3f, d = %.3f",
t.SIFc$parameter, t.SIFc$statistic, t.SIFc$p.value, t.SIFc$statistic/sqrt(t.SIFc$parameter+1))
To identify univariate and bi-variate outliers in the SSRT and SIF scores, we used box plot method, which relies on the interquartile range. Univariate outliers were not present for any of the two measures. One bi-variate outlier was removed from the correlation analysis and the behavioural partial least squares analysis (described below). Nevertheless, outlier removal did not qualitatively alter the results.
# ==================================================================
# SSRT and SIF OUTLIERS
# ==================================================================
SSRT <- dataSSTsubj$SSRT
SIF <- dataTNT$spSIFc
options(repr.plot.width = 7, repr.plot.height = 2, repr.plot.res = 200) # change plot size
u_spipOutliers <- get_outliers(dataTNT$spipSIFu, "SP&IP", "%", SSRT, "SSRT", "ms", disp = TRUE)
u_spOutliers <- get_outliers(dataTNT$spSIFu, "SP", "%", SSRT, "SSRT", "ms", disp = TRUE)
u_ipOutliers <- get_outliers(dataTNT$ipSIFu, "IP", "%", SSRT, "SSRT", "ms", disp = TRUE)
options(repr.plot.width = 7, repr.plot.height = 2, repr.plot.res = 200) # change plot size
c_spipOutliers <- get_outliers(dataTNT$spipSIFc, "SP&IP", "%", SSRT, "SSRT", "ms", disp = TRUE)
c_spOutliers <- get_outliers(dataTNT$spSIFc, "SP", "%", SSRT, "SSRT", "ms", disp = TRUE)
c_ipOutliers <- get_outliers(dataTNT$ipSIFc, "IP", "%", SSRT, "SSRT", "ms", disp = TRUE)
options(repr.plot.width = 11, repr.plot.height = 6, repr.plot.res = 200) # change plot size
grid.arrange(
arrangeGrob(
plotCorrelation(dataTNT$spipSIFu,SSRT,"Inhibiting memories (SP&IP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$spipSIFu, SSRT, u_spipOutliers),
1.8, 11, u_spipOutliers),
plotCorrelation(dataTNT$spSIFu,SSRT,"Inhibiting memories (SP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$spSIFu, SSRT, u_spOutliers),
1.8, 11, u_spOutliers),
plotCorrelation(dataTNT$ipSIFu,SSRT,"Inhibiting memories (IP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$ipSIFu, SSRT, u_ipOutliers),
1.8, 11, u_ipOutliers),
top = textGrob("Unconditionalised", gp=gpar(fontface="bold")),
nrow = 1),
arrangeGrob(
plotCorrelation(dataTNT$spipSIFc,SSRT,"Inhibiting memories (SP&IP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$spipSIFc, SSRT, c_spipOutliers),
1.8, 11, c_spipOutliers) +
theme(panel.background = element_rect(fill = "lightyellow")),
plotCorrelation(dataTNT$spSIFc,SSRT,"Inhibiting memories (SP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$spSIFc, SSRT, c_spOutliers),
1.8, 11, u_spOutliers),
plotCorrelation(dataTNT$ipSIFc,SSRT,"Inhibiting memories (IP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$ipSIFc, SSRT, c_ipOutliers),
1.8, 11, u_ipOutliers),
top = textGrob("Conditionalised", gp=gpar(fontface="bold")),
nrow = 1)
)
plot.behcor <- plotCorrelation(dataTNT$spipSIFc,SSRT,
"Suppression induced forgetting (%)",
"Stop-signal reaction time (ms)",
doCorrelation(dataTNT$spipSIFc, SSRT, c_spipOutliers),
1.8, 11, c_spipOutliers)
# display
options(repr.plot.width = 5, repr.plot.height = 4, repr.plot.res = 250)
plot.behcor
#save
ggsave(file="Figure4a.svg", plot.behcor, width=4, height=3)
# ==================================================================
# SSRT with Trigger Failures
# ==================================================================
# Descriptives
# ----------------------------------------------------------------------
dscr_tfSSRT <-
stat.desc(data_tfSSRT,
basic = TRUE,
desc = TRUE,
p = 0.95)
print(round(dscr_tfSSRT, 2))
tfSSRT <- data_tfSSRT$tfSSRT
tf <- data_tfSSRT$p_tf
options(repr.plot.width = 7, repr.plot.height = 2) # change plot size
u_spipOutl_tfSSRT <- get_outliers(dataTNT$spipSIFu, "SP&IP", "%", tfSSRT, "tfSSRT", "ms", disp = TRUE)
u_spOutl_tfSSRT <- get_outliers(dataTNT$spSIFu, "SP", "%", tfSSRT, "tfSSRT", "ms", disp = TRUE)
u_ipOutl_tfSSRT <- get_outliers(dataTNT$ipSIFu, "IP", "%", tfSSRT, "tfSSRT", "ms", disp = TRUE)
options(repr.plot.width = 7, repr.plot.height = 2) # change plot size
c_spipOutl_tfSSRT <- get_outliers(dataTNT$spipSIFc, "SP&IP", "%", tfSSRT, "tfSSRT", "ms", disp = TRUE)
c_spOutl_tfSSRT <- get_outliers(dataTNT$spSIFc, "SP", "%", tfSSRT, "tfSSRT", "ms", disp = TRUE)
c_ipOutl_tfSSRT <- get_outliers(dataTNT$ipSIFc, "IP", "%", tfSSRT, "tfSSRT", "ms", disp = TRUE)
options(repr.plot.width = 11, repr.plot.height = 6, repr.plot.res = 200) # change plot size
grid.arrange(
arrangeGrob(
plotCorrelation(dataTNT$spipSIFu,tfSSRT,"Inhibiting memories (SP&IP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$spipSIFu, tfSSRT, u_spipOutl_tfSSRT),
1.8, 11, u_spipOutl_tfSSRT),
plotCorrelation(dataTNT$spSIFu,tfSSRT,"Inhibiting memories (SP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$spSIFu, tfSSRT, u_spOutl_tfSSRT),
1.8, 11, u_spOutl_tfSSRT),
plotCorrelation(dataTNT$ipSIFu,tfSSRT,"Inhibiting memories (IP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$ipSIFu, tfSSRT, u_ipOutl_tfSSRT),
1.8, 11, u_ipOutl_tfSSRT),
top = textGrob("Unconditionalised", gp=gpar(fontface="bold")),
nrow = 1),
arrangeGrob(
plotCorrelation(dataTNT$spipSIFc,tfSSRT,"Inhibiting memories (SP&IP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$spipSIFc, tfSSRT, c_spipOutl_tfSSRT),
1.8, 11, c_spipOutl_tfSSRT) +
theme(panel.background = element_rect(fill = "lightyellow")),
plotCorrelation(dataTNT$spSIFc,tfSSRT,"Inhibiting memories (SP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$spSIFc, tfSSRT, c_spOutl_tfSSRT),
1.8, 11, u_spOutl_tfSSRT),
plotCorrelation(dataTNT$ipSIFc,tfSSRT,"Inhibiting memories (IP SIF, %)",
"Inhibiting actions (SSRT, ms)",
doCorrelation(dataTNT$ipSIFc, tfSSRT, c_ipOutl_tfSSRT),
1.8, 11, u_ipOutl_tfSSRT),
top = textGrob("Conditionalised", gp=gpar(fontface="bold")),
nrow = 1)
)
options(repr.plot.width = 5, repr.plot.height = 3, repr.plot.res = 200) # change plot size
ssrt.tf.outl <- get_outliers(SSRT, "SSRT", "%", tf, "TF", "p", disp = FALSE)
plotCorrelation(SSRT,tf,"SSRT","Stop trigger failures (p)",
doCorrelation(SSRT, tf, ssrt.tf.outl),
1.8, 11, ssrt.tf.outl)
options(repr.plot.width = 7, repr.plot.height = 2, repr.plot.res = 200) # change plot size
u_spipOutl_tf <- get_outliers(dataTNT$spipSIFu, "SP&IP", "%", tf, "TF", "p", disp = FALSE)
u_spOutl_tf <- get_outliers(dataTNT$spSIFu, "SP", "%", tf, "TF", "p", disp = FALSE)
u_ipOutl_tf <- get_outliers(dataTNT$ipSIFu, "IP", "%", tf, "TF", "p", disp = FALSE)
options(repr.plot.width = 7, repr.plot.height = 2, repr.plot.res = 200) # change plot size
c_spipOutl_tf <- get_outliers(dataTNT$spipSIFc, "SP&IP", "%", tf, "TF", "p", disp = FALSE)
c_spOutl_tf <- get_outliers(dataTNT$spSIFc, "SP", "%", tf, "TF", "p", disp = FALSE)
c_ipOutl_tf <- get_outliers(dataTNT$ipSIFc, "IP", "%", tf, "TF", "p", disp = FALSE)
options(repr.plot.width = 11, repr.plot.height = 6, repr.plot.res = 200) # change plot size
grid.arrange(
arrangeGrob(
plotCorrelation(dataTNT$spipSIFu,tf,"Inhibiting memories (SP&IP SIF, %)","Stop trigger failures (p)",
doCorrelation(dataTNT$spipSIFu, tf, u_spipOutl_tf),
1.8, 11, u_spipOutl_tf),
plotCorrelation(dataTNT$spSIFu,tf,"Inhibiting memories (SP SIF, %)","Stop trigger failures (p)",
doCorrelation(dataTNT$spSIFu, tf, u_spOutl_tf),
1.8, 11, u_spOutl_tf),
plotCorrelation(dataTNT$ipSIFu,tf,"Inhibiting memories (IP SIF, %)","Stop trigger failures (p)",
doCorrelation(dataTNT$ipSIFu, tf, u_ipOutl_tf),
1.8, 11, u_ipOutl_tf),
top = textGrob("Unconditionalised", gp=gpar(fontface="bold")),
nrow = 1),
arrangeGrob(
plotCorrelation(dataTNT$spipSIFc,tf,"Inhibiting memories (SP&IP SIF, %)","Stop trigger failures (p)",
doCorrelation(dataTNT$spipSIFc, tf, c_spipOutl_tf),
1.8, 11, c_spipOutl_tf) + theme(panel.background = element_rect(fill = "lightyellow")),
plotCorrelation(dataTNT$spSIFc,tf,"Inhibiting memories (SP SIF, %)","Stop trigger failures (p)",
doCorrelation(dataTNT$spSIFc, tf, c_spOutl_tf),
1.8, 11, u_spOutl_tf),
plotCorrelation(dataTNT$ipSIFc,tf,"Inhibiting memories (IP SIF, %)","Stop trigger failures (p)",
doCorrelation(dataTNT$ipSIFc, tf, c_ipOutl_tf),
1.8, 11, u_ipOutl_tf),
top = textGrob("Conditionalised", gp=gpar(fontface="bold")),
nrow = 1)
)
# ----------------------------------------------------------------------
# GET THE DATA
# ----------------------------------------------------------------------
# Load
plsdata <-
read.csv('https://raw.githubusercontent.com/dcdace/tmp/master/data/pls_brain_scores.csv')
plsdata <- plsdata
plsdata$SSRT <- SSRT
plsdata$SIF <- dataTNT$spipSIFc
pls.ssrt.outl <- get_outliers(plsdata$SSRT, "SSRT", "ms",
plsdata$BrainScores, "Brain scores", "", disp = FALSE)
pls.sif.outl <- get_outliers(plsdata$SIF, "SIF", "%",
plsdata$BrainScores, "Brain scores", "", disp = FALSE)
plot.pls.ssrt <- plotCorrelation(plsdata$SSRT, plsdata$BrainScores,
"SSRT", "Brain scores",
doCorrelation(plsdata$SSRT, plsdata$BrainScores, pls.ssrt.outl),
1.8, 14, pls.ssrt.outl, FALSE)
plot.pls.sif <- plotCorrelation(plsdata$SIF, plsdata$BrainScores,
"SIF", "Brain scores",
doCorrelation(plsdata$SIF, plsdata$BrainScores, pls.sif.outl),
1.8, 14, pls.sif.outl, FALSE)
options(repr.plot.width = 9, repr.plot.height = 3, repr.plot.res = 200)
grid.arrange(plot.pls.ssrt, plot.pls.sif, nrow = 1)
# ----------------------------------------------------------------------
# GET THE DATA
# ----------------------------------------------------------------------
# Load
data <-
read.csv('https://raw.githubusercontent.com/dcdace/tmp/master/data/targetROIs_psc.csv')
data[, 4:9] = data[, 4:9]-data[, 4]
# Have a look what it contains
head(data,1)
cat('rois: ')
cat(unique(data$roi))
cat('\nconditions: ')
cat(unique(data$condition))
# ----------------------------------------------------------------------
# ROI x MODALITY INTERACTION
# ----------------------------------------------------------------------
# Create Inhibit (NT & Stop) and Express (T & Go) subsets
dataInh <- subset(data, (condition == 'nt' | condition == 's'))
dataExp <- subset(data, (condition == 't' | condition == 'g'))
# Substract Express from Inhibit
dataModality <- cbind(dataInh[, 1:3], dataInh[, 4:9] - dataExp[, 4:9])
dataModality$condition <- as.factor(dataModality$condition)
# Change conditions to the corresponding Modality
dataModality$condition <-
revalue(dataModality$condition, c("nt" = "nt_t", "s" = "s_g"))
names(dataModality)[names(dataModality) == 'condition'] <-
'modality'
# Average PSC within a timespan of 2-6s and substracting the onset value to account for pretrial variability.
# Similar as in Levy&Anderson, 2012 (they did 4-8s)
dataModality$PSC <- rowMeans(dataModality[, 5:7]) #-dataModality[, 4]
# Change the Type levels to better reflect the data content
dataModality$modality <- factor(dataModality$modality,
levels = c("nt_t", "s_g"),
labels = c("No-Think - Think", "Stop - Go"))
# See what dataModality contains now
head(dataModality,1)
levels(dataModality$modality)
# Using rm_2by2_anova.R function from https://raw.githubusercontent.com/dcdace/R_functions/
# Specify columns of interest
columns <- list(
sID = "sid", # subject ID
DV = "PSC", # dependent variable
Fc1 = "modality", # within-subj factor 1
Fc2 = "roi" # within-subj factor 2
)
# Define plot label names, tilte, colors and error bar type
param <- list(
y.label = "% signal change difference",
Fc1.label = "Modality",
Fc2.label = "Target region",
cat.color = c('red', '#ffc000'),
errorbar = 'se'
)
# the plot title
param$title <- sprintf('%s x %s interaction', param$Fc2.label, param$Fc1.label)
# Run the ANOVA function but don't display the printouts
anova.results <- hush(rm_2by2_anova(dataModality, columns, param))
# Adding post-hoc comparison stars
# Prepare the values
# post-hoc pairwise comparisons for each Fc1 level between the Fc2 levels
p.Fc1.L1 <-
ifelse(anova.results$pwc1$p[1] < 0.05, stars.pval(anova.results$pwc1$p[1]), "n.s.")
p.Fc1.L2 <-
ifelse(anova.results$pwc1$p[2] < 0.05, stars.pval(anova.results$pwc1$p[2]), "n.s.")
# add horizontal line behind the plot
anova.results$plot.anova$layers <- c(geom_hline(aes(yintercept = 0), size = 0.1), # the new layer
anova.results$plot.anova$layers) # original layers
# add other things
plot.targetrois <- anova.results$plot.anova +
# Add post-hoc pairwise comparisons
# might need to adjust the x and y positions
annotate(
"text",
x = 1.38,
y = -0.06,
label = p.Fc1.L1,
color = param$cat.color[1],
size = 8
) +
annotate(
"text",
x = 1.6,
y = -0.15,
label = p.Fc1.L2,
color = param$cat.color[2],
size = 8
) +
# Change ROI captions
scale_x_discrete(limits = c("rHC", "lM1"), label = c('Hippocampus', 'M1'))
#display
options(repr.plot.width = 5, repr.plot.height = 4, repr.plot.res = 200) # change plot size
plot.targetrois
#save
ggsave(file="Figure5a.pdf", plot.targetrois, width=5, height=4)
summary.anova <- summary(anova.results$res.anova)
res.txt.anova <- sprintf(
'%s: F(%d,%d) = %.2f, p = %.7f',
param$title,
summary.anova[["Error: sID:Fc1:Fc2"]][[1]][["Df"]][1],
summary.anova[["Error: sID:Fc1:Fc2"]][[1]][["Df"]][2],
summary.anova[["Error: sID:Fc1:Fc2"]][[1]][["F value"]][1],
summary.anova[["Error: sID:Fc1:Fc2"]][[1]][["Pr(>F)"]][1]
)
cat('\n================================================\n\n')
cat(res.txt.anova)
cat('\n\n================================================')
# Main effect of Modality
cat('Main effect of Modality')
anova.results$one.way1
# Main effect of Target region
cat('Main effect of Target region')
anova.results$one.way2
res.txt.ntt <- sprintf(
't(%d) = %.2f, p = %.3f, d = %.3f',
anova.results$pwc1[1,]$df,
abs(anova.results$pwc1[1,]$statistic),
anova.results$pwc1[1,]$p,
abs(anova.results$pwc1[1,]$statistic / sqrt(anova.results$pwc1[1,]$df + 1))
)
cat(res.txt.ntt)
res.txt.sg <- sprintf(
't(%d) = %.2f, p = %.5f, d = %.3f',
anova.results$pwc1[2,]$df,
abs(anova.results$pwc1[2,]$statistic),
anova.results$pwc1[2,]$p,
abs(anova.results$pwc1[2,]$statistic / sqrt(anova.results$pwc1[1,]$df + 1))
)
cat(res.txt.sg)
# Convert to long form
dataLong <- gather(data, seconds, psc, X0s:X10s, factor_key = TRUE)
dataLong$inhibit <-
as.factor(ifelse((
dataLong$condition == 's' | dataLong$condition == 'nt'
), 1, 0))
# Subsets for each of the 4 levels
# rHC T & NT
dataLong.rHCTNT <-
subset(dataLong, roi == 'rHC' &
(condition == 't' | condition == 'nt'))
# rHC Stop & Go
dataLong.rHCSG <-
subset(dataLong, roi == 'rHC' &
(condition == 'g' | condition == 's'))
# lM1 T & NT
dataLong.lM1TNT <-
subset(dataLong, roi == 'lM1' &
(condition == 't' | condition == 'nt'))
# lM1 Stop & Go
dataLong.lM1SG <-
subset(dataLong, roi == 'lM1' &
(condition == 'g' | condition == 's'))
# ----------------------------------------------------------------------
# PLOT FIRs
# ----------------------------------------------------------------------
plotFirs <- function(dataset, labelstxt, title) {
# aggregate
res <-
summarySEwithin(
dataset,
measurevar = "psc",
withinvars = c("seconds", "inhibit"),
idvar = "sid",
na.rm = TRUE,
conf.interval = 0.95
)
# plot
ggplot(data = res, aes(x = seconds, y = psc, group = inhibit)) +
geom_hline(aes(yintercept = 0), size = 0.2) +
geom_line(size = 1.5, aes(color = inhibit)) +
geom_ribbon(
aes(
ymin = psc - se,
ymax = psc + se,
x = seconds,
fill = inhibit
),
alpha = 0.3,
show.legend = FALSE
) +
geom_line(aes(x = seconds, y = psc - se, color = inhibit), size = 0.2) +
geom_line(aes(x = seconds, y = psc + se, color = inhibit), size = 0.2) +
scale_color_manual(labels = labelstxt, values = c('green', 'blue')) +
scale_fill_manual(values = c('#14b814', '#2e599e')) +
scale_x_discrete(labels = c("0", "2", "4", "6", "8", "10")) +
labs(x = "Peristimulus time (s)", y = "% signal change") +
ggtitle(title) +
guides(color = guide_legend("")) +
theme_minimal() +
theme(
legend.position = "top",
text = element_text(size = 20),
plot.title = element_text(
hjust = 0.5,
size = 20,
face = "bold"
),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
}
# ----------------------------------------------------------------------
# M1: Stop/Go Plot
# ----------------------------------------------------------------------
p.gs <- plotFirs(dataLong.lM1SG, c("Go", "Stop"), "Go/Stop, M1")
# Add one-tailed significance for > or < baseline
lm1s <- subset(data, roi == 'lM1' & condition == 's')
lm1s$PSC <- rowMeans(lm1s[, 5:8])
p.lm1s <- stars.pval(sapply(c(5:8), function(x) {
t.test(lm1s[, x], mu = 0, alternative = c("less"))$p.value
}))
lm1g <- subset(data, roi == 'lM1' & condition == 'g')
lm1g$PSC <- rowMeans(lm1s[, 5:8])
p.lm1g <- stars.pval(sapply(c(5:8), function(x) {
t.test(lm1g[, x], mu = 0, alternative = c("greater"))$p.value
}))
# position the significance stars
sstypos <-
summarySEwithin(
dataLong.lM1SG,
measurevar = "psc",
withinvars = c("seconds", "inhibit"),
idvar = "sid",
na.rm = TRUE,
conf.interval = 0.95
)
# add significance stars to the plot
p.gs <- p.gs +
annotate("text", y = min(sstypos$psc) - mean(sstypos$ci),
x = 2:5, label = p.lm1s, size = 10, color = "blue"
) +
annotate(
"text", y = max(sstypos$psc) + mean(sstypos$ci),
x = 2:5, label = p.lm1g, size = 10, color = "green")
# ----------------------------------------------------------------------
# Hippocampus: T/NT Plot
# ----------------------------------------------------------------------
p.tnt <- plotFirs(dataLong.rHCTNT, c("Think", "No-Think"), "Think/No-Think, Hippocampus")
# Add one-tailed significance for > or < baseline
rhcnt <- subset(data, roi == 'rHC' & condition == 'nt')
rhcnt$PSC <- rowMeans(rhcnt[, 5:8])
p.rhcnt <- stars.pval(sapply(c(5:8), function(x) {
t.test(rhcnt[, x], mu = 0, alternative = c("less"))$p.value
}))
rhct <- subset(data, roi == 'rHC' & condition == 't')
rhct$PSC <- rowMeans(rhct[, 5:8])
p.rhct <- stars.pval(sapply(c(5:8), function(x) {
t.test(rhct[, x], mu = 0, alternative = c("greater"))$p.value
}))
# position the significance stars
tntypos <-
summarySEwithin(
dataLong.rHCTNT,
measurevar = "psc",
withinvars = c("seconds", "inhibit"),
idvar = "sid",
na.rm = TRUE,
conf.interval = 0.95)
# add significance stars to the plot
p.tnt <- p.tnt +
annotate(
"text", y = min(tntypos$psc) - mean(tntypos$ci),
x = 2:5, label = p.rhcnt,size = 10, color = "blue"
) +
annotate("text", y = max(tntypos$psc) + mean(tntypos$ci),
x = 2:5, label = p.rhct, size = 10, color = "green"
)
#display
options(repr.plot.width = 15, repr.plot.height = 5, repr.plot.res = 200)
plot.firs <- grid.arrange(p.gs, p.tnt, nrow=1)
#save
ggsave(file="Figure5b_gs.pdf", p.gs, width = 6, height = 3)
ggsave(file="Figure5b_tnt.pdf", p.tnt, width = 6, height = 3)
data.s.lM1 = subset(data, condition == "s" & roi == "lM1")
data.s.lM1$PSC <- rowMeans(data.s.lM1[, 5:7]) - data.s.lM1[, 4]
tt_sm1 <- t.test(data.s.lM1$PSC, mu = 0, alternative = "less")
cat(sprintf('t(%d) = %.2f, p = %.4f, d = %.3f',
tt_sm1$parameter,
tt_sm1$statistic,
tt_sm1$p.value,
abs(tt_sm1$statistic/sqrt(tt_sm1$parameter+1))
))
data.nt.rHC = subset(data, condition == "nt" & roi == "rHC")
data.nt.rHC$PSC <- rowMeans(data.nt.rHC[, 5:7]) - data.nt.rHC[, 4]
tt_nthc <- t.test(data.nt.rHC$PSC, mu = 0, alternative = "less")
cat(sprintf('t(%d) = %.2f, p = %.3f, d = %.3f',
tt_nthc$parameter,
tt_nthc$statistic,
tt_nthc$p.value,
abs(tt_nthc$statistic/sqrt(tt_nthc$parameter+1))
))
# ----------------------------------------------------------------------
# GET THE DATA
# ----------------------------------------------------------------------
df <-
read.csv(
'https://raw.githubusercontent.com/dcdace/tmp/master/data/cross_decoding.csv',
fileEncoding = "UTF-8-BOM"
)
df$acc <- df$acc*100
df$roi <- factor(df$roi, levels=c('rdlpfc', 'rvlpfc', 'rHC', 'lM1'))
#df$type <- factor(df$type, levels=c('d_specific', 'd_general'))
rois <- levels(df$roi)
nrois <- length(rois)
# CREATE SUBSETS
# ------------------------------------
# Action domain decoding
action <- subset(df, df$type == 'ss')
# Thought domain decoding
thought <- subset(df, df$type == 'ntnt')
# Within-domain decoding
# 1) Train: NT, T; Test: NT=NT
# 2) Train: Stop, Go; Test: Stop=Stop
# average both
within <- subset(df, df$type == 'within_d')
# Inhibition decoding, domain-special
# Train: NT, Stop; Test: NT<->Stop
specific <- subset(df, df$type == 'd_specific')
# Cross-decoding, domain-general
# 1) Train: Stop,Go; Test: NT,T
# 2) Train: NT,T; Test: Stop,Go
# average both
general <- subset(df, df$type == 'd_general')
# Cross Action-to-Thought
# Train: Stop,Go; Test: NT=Stop
a_to_t <- subset(df, df$type == 'action_to_thought')
# Cross Thought-to-Action
# Train: NT,T; Test: Stop=NT
t_to_a <- subset(df, df$type == 'thought_to_action')
# TWO DOMAIN PLOTS AND STATS
two_domain_analysis <-
function(data, color1, color2, name1, name2, plottitle) {
# Within-Subject summary for errorbars
dataSummary <- summarySEwithin(
data,
measurevar = "acc",
withinvars = c("roi", "type"),
idvar = "sid",
na.rm = TRUE,
conf.interval = 0.95
)
# PLOT
# ------------------------------------
plot.mvpa <-
ggplot(data = data, aes(x = roi, y = acc, group = type)) +
geom_hline(aes(yintercept = 50), size = 0.2) +
geom_point(
colour = "black", alpha = .5, aes(fill = type), size = 2, stroke = 0.5, shape = 21,
position = position_jitterdodge()
) +
geom_pointrange(
data = dataSummary,
aes(y = acc, ymin = acc - se, ymax = acc + se, fill = type
),
alpha = 1, color = "black", size = 1.5, stroke = 0.5, shape = 21,
position = position_dodge(width = .2),
show.legend = FALSE
) +
scale_fill_manual(
values = c(color1, color2),
name = "",
labels = c(name1, name2),
guide = guide_legend(
override.aes = list(
alpha = 1, color = "black", size = 5.5, stroke = 0.5, shape = 21
)
)
) +
scale_y_continuous(limits = c(min(data$acc) - sd(data$acc), max(data$acc)),
breaks = seq(0, 100, 20)) +
labs(x = "", y = "Classification accuracy (%)") + ggtitle(plottitle) +
scale_x_discrete(labels = c("rDLPFC", "rVLPFC", "Hippocampus", "M1")) +
theme_minimal() +
theme(
text = element_text(size = 22),
plot.title = element_text(
size = 20, face = "bold", hjust = 0.5
),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "top"
)
# ADD STATS
# ------------------------------------
# One-sample t-test, > 50%, Bonf. corrected for 4 ROIs
stat.test <- data %>%
group_by(roi, type) %>%
t_test(acc ~ 1, mu = 50, alternative = "greater") %>%
adjust_pvalue() %>%
mutate(y.position = 20)
stat.test$stars <-
stars.pval(ifelse(stat.test$p * 4 > 1, 1, stat.test$p * 4))
stat.test$colors <-
ifelse(stat.test$type == "d_general" |
stat.test$type == "action_to_thought" |
stat.test$type == "thought_to_action",
color2, color1)
plot.mvpa <- plot.mvpa +
stat_pvalue_manual(
stat.test,
label = "stars",
x = "roi",
position = position_dodge(0.8),
size = 10,
color = stat.test$colors
)
return(plot.mvpa)
}
# ------------------------------------
# Paired t-tests, Bonf. corrected for 4 ROIs
# paired ttest of domain-general vs domain-specific
pairedttest <- function(data1, data2){
test <- t.test(data1$acc, data2$acc, paired = TRUE)
# correct p-value for number of rois tested
test$p.value <- ifelse(test$p.value * nrois > 1, 1, test$p.value * nrois)
return(
data.frame(
tval = test$statistic,
df = test$parameter[1],
pval = test$p.value,
M = mean(data1$acc) - mean(data2$acc)
))
}
# ------------------------------------
# t-test function
dottest <- function(data) {
test <- t.test(data$acc, mu = 50, alternative = 'greater')
# correct p-value for number of rois tested and keep max p to 1
test$p.value <-
ifelse(test$p.value * nrois > 1, 1, test$p.value * nrois)
return(data.frame(
tval = test$statistic,
df = test$parameter[1],
pval = test$p.value,
M = mean(data$acc),
std = sd(data$acc)
))
}
dispresults <- function(res) {
txt <- lapply(seq(1:nrois), function(r) {
sprintf(
"%s: M = %.0f%% (+-%.0f), t(%d) = %.2f, p = %.3f, d = %.3f",
res$roi[r],
res$M[r],
res$std[r],
res$df[r],
res$tval[r],
res$pval[r],
res$tval[r] / sqrt(res$df[r] + 1)
)
})
txtBind <- (ldply(txt, "rbind", .id = "roi"))
}
# ------------------------------------
# PLOT WITHIN-SUBJECT MEANS
# ------------------------------------
plot_within_subj_means <-
function(data, yvar, groupvar, sid, ylabel, xlabel,title) {
# data summary to get within SE and CI
dataSummary <- summarySEwithin(
data,
measurevar = yvar,
withinvars = groupvar,
idvar = sid,
na.rm = TRUE,
conf.interval = 0.95
)
# plot
ggplot(data, aes(data$roi, data$acc)) +
geom_point(
colour = "black", alpha = .9, fill = 'grey', size = 1, stroke = 0.1, shape = 21,
position = position_jitter(height = 0, width = 0.1)
) +
geom_pointrange(
data = dataSummary,
aes(
y = dataSummary$acc,
ymin = dataSummary$acc - se,
ymax = dataSummary$acc + se,
x = dataSummary$roi
),
fill = 'red', alpha = 0.8, color = "black", size = 0.6, stroke = 0.3,
shape = 21,position = position_nudge(x = -0.2)
) +
labs(x = xlabel, y = ylabel) +
ggtitle(title) +
scale_x_discrete(label = c('rDLPFC', 'rVLPFC', 'Hippocampus', 'M1')) +
geom_hline(aes(yintercept = 0), size = 0.1) +
theme_minimal() +
theme(
text = element_text(size = 9),
plot.title = element_text(
hjust = 0.5, size = 9, face = "bold"
),
plot.background = element_rect(
color = "lightgrey", size = 0.4, linetype = "dotted"
),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
}
Domain-specific component represents how well the classifier discriminates between Stop and No-Think conditions, Stop $\ne$ Think
Domain-general component represents how well the classifier discriminates between Inhibition and Non-Inhibition conditions independent of modality (training the classifier to discriminate No-Think from Think and testing how well it discriminates Stop from Go, and vice versa). Cross-modality Inhibit $\ne$ Non-Inhibit.
#
df_dd <- subset(df, type == 'd_general' | type == 'd_specific')
df_dd$type <- factor(df_dd$type, levels = c('d_specific', 'd_general'))
# PLOT
# ------------------------------------
plot_dd <- two_domain_analysis(df_dd,
"#aad400ff", '#ff00ff',
"Domain-specific component", "Domain-general component",
"Inhibition representations")
# do ttests per ROI for each subset of dataset
# ------------------------------------
tres <- sapply(list(specific, general), function(x) {
ddply(x, .(roi), dottest)
})
resdf <-
as.data.frame(
cbind(
dispresults(tres[, 1]),
dispresults(tres[, 2])
)
)
colnames(resdf) <-
c("domain-specific", "domain-general")
resdf
# Paired-ttest
# ------------------------------------
ptests <- sapply(seq(1:nrois), function(x) {
pairedttest(subset(specific, roi == rois[x]), subset(general, roi == rois[x]))
})
# get stars for pvalues
sig <- sapply(seq(1:nrois), function(x) {
stars.pval(ptests[, x]$pval)
})
#add significance stars to the plot
plot_dd <- plot_dd +
annotate(
"text",
y = 18,
x = 1:4,
label = sig,
size = 10,
color = "darkgrey"
)
txt <- lapply(seq(1:nrois), function(x) {
sprintf(
"%s: M = %.0f%%, t(%d) = %.2f, p = %.3f, d = %.3f",
rois[x],
ptests[, x]$M,
ptests[, x]$df,
ptests[, x]$tval,
ptests[, x]$pval,
ptests[, x]$tval / sqrt(ptests[, x]$df + 1)
)
})
ldply(txt, "rbind", .id = "roi")
options(repr.plot.width = 10, repr.plot.height = 8, repr.plot.res = 150)
plot_dd
#save
ggsave(file="Figure5c.pdf", plot_dd, width = 10, height = 8)
# DIFFERENCE
# ------------------------------------
# prepare data
sg_dif <- specific[, c("sid", "roi")]
sg_dif$acc <- specific$acc - general$acc
# 1-way RM ANOVA
res.aov <-
anova_test(data = sg_dif, dv = acc, wid = sid, within = roi)
get_anova_table(res.aov)
# pairwise comparisons
pwc <- sg_dif %>%
pairwise_t_test(acc ~ roi, paired = TRUE,
p.adjust.method = "bonferroni")
pwc
# plot
cat("Error bars: SE")
options(repr.plot.width = 3, repr.plot.height = 3, repr.plot.res = 300)
p_dd_dif <- plot_within_subj_means(sg_dif, "acc", "roi", "sid",
"Classification accuracy difference", "",
"Domain-specific > Domain-general")
p_dd_dif
#save
ggsave(file="FigureR1.3.dd.pdf", p_dd_dif, width = 3, height = 3)
Within-domain: an average of Stop $\ne$ Go and No-Think $\ne$ Think
#
df_w <- subset(df, type == 'd_general' | type == 'within_d')
df_w$type <- factor(df_w$type, levels = c('within_d', 'd_general'))
# PLOT
# ------------------------------------
plot_w <- two_domain_analysis(df_w,
"#468ce2", '#ff00ff',
"Within-domain", "Domain-general",
"Inhibition representations")
# do ttests per ROI for each subset of dataset
# ------------------------------------
tres <- sapply(list(within, general), function(x) {
ddply(x, .(roi), dottest)
})
resdf <-
as.data.frame(
cbind(
dispresults(tres[, 1]),
dispresults(tres[, 2])
)
)
colnames(resdf) <-
c("within-domain", "domain-general")
resdf
# Paired-ttest
# ------------------------------------
ptests <- sapply(seq(1:nrois), function(x) {
pairedttest(subset(within, roi == rois[x]), subset(general, roi == rois[x]))
})
# get stars for pvalues
sig <- sapply(seq(1:nrois), function(x) {
stars.pval(ptests[, x]$pval)
})
#add significance stars to the plot
plot_w <- plot_w +
annotate(
"text",
y = 18,
x = 1:4,
label = sig,
size = 10,
color = "darkgrey"
)
txt <- lapply(seq(1:nrois), function(x) {
sprintf(
"%s: M = %.0f%%, t(%d) = %.2f, p = %.3f, d = %.3f",
rois[x],
ptests[, x]$M,
ptests[, x]$df,
ptests[, x]$tval,
ptests[, x]$pval,
ptests[, x]$tval / sqrt(ptests[, x]$df + 1)
)
})
ldply(txt, "rbind", .id = "roi")
options(repr.plot.width = 10, repr.plot.height = 8, repr.plot.res = 150)
plot_w
#save
ggsave(file="Figure5c_within.pdf", plot_w, width = 10, height = 8)
# DIFFERENCE
# ------------------------------------
# prepare data
sg_dif <- within[, c("sid", "roi")]
sg_dif$acc <- within$acc - general$acc
# 1-way RM ANOVA
res.aov <-
anova_test(data = sg_dif, dv = acc, wid = sid, within = roi)
get_anova_table(res.aov)
# pairwise comparisons
pwc <- sg_dif %>%
pairwise_t_test(acc ~ roi, paired = TRUE,
p.adjust.method = "bonferroni")
pwc
# plot
cat("Error bars: SE")
options(repr.plot.width = 3, repr.plot.height = 3, repr.plot.res = 300)
p_w_dif <- plot_within_subj_means(sg_dif, "acc", "roi", "sid",
"Classification accuracy difference", "",
"Within-domain > Domain-general")
p_w_dif
#save
ggsave(file="FigureR1.3.w.pdf", p_w_dif, width = 3, height = 3)
Action domain: discriminates between Stop and Go, Stop $\ne$ Go
#
df_a <- subset(df, type == 'action_to_thought' | type == 'ss')
df_a$type <- factor(df_a$type, levels = c('ss', 'action_to_thought'))
# PLOT
# ------------------------------------
plot_a <- two_domain_analysis(df_a,
"#ffc000", '#ff00ff',
"Action-domain", "Cross: Action-to-Thought",
"Inhibition representations")
# do ttests per ROI for each subset of dataset
# ------------------------------------
tres <- sapply(list(action, a_to_t), function(x) {
ddply(x, .(roi), dottest)
})
resdf <-
as.data.frame(
cbind(
dispresults(tres[, 1]),
dispresults(tres[, 2])
)
)
colnames(resdf) <-
c("action-domain", "cross: action-to-thought")
resdf
# Paired-ttest
# ------------------------------------
ptests <- sapply(seq(1:nrois), function(x) {
pairedttest(subset(action, roi == rois[x]), subset(a_to_t, roi == rois[x]))
})
# get stars for pvalues
sig <- sapply(seq(1:nrois), function(x) {
stars.pval(ptests[, x]$pval)
})
#add significance stars to the plot
plot_a <- plot_a +
annotate(
"text",
y = 18,
x = 1:4,
label = sig,
size = 10,
color = "darkgrey"
)
txt <- lapply(seq(1:nrois), function(x) {
sprintf(
"%s: M = %.0f%%, t(%d) = %.2f, p = %.3f, d = %.3f",
rois[x],
ptests[, x]$M,
ptests[, x]$df,
ptests[, x]$tval,
ptests[, x]$pval,
ptests[, x]$tval / sqrt(ptests[, x]$df + 1)
)
})
ldply(txt, "rbind", .id = "roi")
options(repr.plot.width = 10, repr.plot.height = 8, repr.plot.res = 150)
plot_a
#save
ggsave(file="Figure5c_action.pdf", plot_a, width = 10, height = 8)
# DIFFERENCE
# ------------------------------------
# prepare data
sg_dif <- action[, c("sid", "roi")]
sg_dif$acc <- action$acc - a_to_t$acc
# 1-way RM ANOVA
res.aov <-
anova_test(data = sg_dif, dv = acc, wid = sid, within = roi)
get_anova_table(res.aov)
# pairwise comparisons
pwc <- sg_dif %>%
pairwise_t_test(acc ~ roi, paired = TRUE,
p.adjust.method = "bonferroni")
pwc
# plot
cat("Error bars: SE")
options(repr.plot.width = 3, repr.plot.height = 3, repr.plot.res = 300)
p_a_dif <- plot_within_subj_means(sg_dif, "acc", "roi", "sid",
"Classification accuracy difference", "",
"Action-specific > Cross: Action-to-Thought")
p_a_dif
#save
ggsave(file="FigureR1.3.a.pdf", p_a_dif, width = 3, height = 3)
Thought domain: discriminates between No-Think and Think, No-Think $\ne$ Think
#
df_t <- subset(df, type == 'thought_to_action' | type == 'ntnt')
df_t$type <- factor(df_t$type, levels = c('ntnt', 'thought_to_action'))
# PLOT
# ------------------------------------
plot_t <- two_domain_analysis(df_t,
"red", '#ff00ff',
"Thought-domain", "Cross: Thought-to-Action",
"Inhibition representations")
# do ttests per ROI for each subset of dataset
# ------------------------------------
tres <- sapply(list(thought, t_to_a), function(x) {
ddply(x, .(roi), dottest)
})
resdf <-
as.data.frame(
cbind(
dispresults(tres[, 1]),
dispresults(tres[, 2])
)
)
colnames(resdf) <-
c("thought-domain", "cross: thought_to_action")
resdf
# Paired-ttest
# ------------------------------------
ptests <- sapply(seq(1:nrois), function(x) {
pairedttest(subset(thought, roi == rois[x]), subset(t_to_a, roi == rois[x]))
})
# get stars for pvalues
sig <- sapply(seq(1:nrois), function(x) {
stars.pval(ptests[, x]$pval)
})
#add significance stars to the plot
plot_t <- plot_t +
annotate(
"text",
y = 18,
x = 1:4,
label = sig,
size = 10,
color = "darkgrey"
)
txt <- lapply(seq(1:nrois), function(x) {
sprintf(
"%s: M = %.0f%%, t(%d) = %.2f, p = %.3f, d = %.3f",
rois[x],
ptests[, x]$M,
ptests[, x]$df,
ptests[, x]$tval,
ptests[, x]$pval,
ptests[, x]$tval / sqrt(ptests[, x]$df + 1)
)
})
ldply(txt, "rbind", .id = "roi")
options(repr.plot.width = 10, repr.plot.height = 8, repr.plot.res = 150)
plot_t
#save
ggsave(file="Figure5c_thought.pdf", plot_t, width = 10, height = 8)
# DIFFERENCE
# ------------------------------------
# prepare data
sg_dif <- thought[, c("sid", "roi")]
sg_dif$acc <- thought$acc - t_to_a$acc
# 1-way RM ANOVA
res.aov <-
anova_test(data = sg_dif, dv = acc, wid = sid, within = roi)
get_anova_table(res.aov)
# pairwise comparisons
pwc <- sg_dif %>%
pairwise_t_test(acc ~ roi, paired = TRUE,
p.adjust.method = "bonferroni")
pwc
# plot
cat("Error bars: SE")
options(repr.plot.width = 3, repr.plot.height = 3, repr.plot.res = 300)
p_t_dif <- plot_within_subj_means(sg_dif, "acc", "roi", "sid",
"Classification accuracy difference", "",
"Domain-specific > Cross: Thought-to-Action")
p_t_dif
#save
ggsave(file="FigureR1.3.t.pdf", p_t_dif, width = 3, height = 3)
# ----------------------------------------------------------------------
# PLOTTING FUNCTION
# ----------------------------------------------------------------------
plotAcc <- function(dataset, labelstxt, title) {
# aggregate
res1 <- dataset %>% # this is to get correct means, ignoring NaNs
group_by(run) %>%
get_summary_stats(acc, type = "mean")
res2 <- # this is to get within-subject SEs
summarySEwithin(
dataset,
measurevar = "acc",
withinvars = c("run"),
idvar = "sNR",
na.rm = TRUE,
conf.interval = 0.95,
)
res <- data.frame(run = res1$run, acc = res1$mean, se = res2$se)
# plot
ggplot(data = res, aes(x = run, y = acc, group=1)) +
geom_hline(aes(yintercept = 50), size = 0.2, color = "red") +
geom_line(size = 0.5) +
geom_point() +
geom_ribbon(
aes(
ymin = acc - se,
ymax = acc + se,
x = run
),
alpha = 0.3,
fill = 'grey',
show.legend = FALSE
) +
ylim(20,100) +
#geom_line(aes(x = run, y = acc - se), size = 0.1) +
#geom_line(aes(x = run, y = acc + se), size = 0.1) +
labs(x = "Run", y = "Classification accuracy (%)") +
geom_text(aes(label = sprintf('%d%%', round(acc,0))), vjust = -1.5) +
ggtitle(title) +
guides(color = guide_legend("")) +
theme_minimal() +
theme(
legend.position = "top",
text = element_text(size = 14),
plot.title = element_text(
hjust = 0.5,
size = 14,
face = "bold"
),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
}
## Per-sun accuracies
perRunResults <- function(data){
Slope <- data$propslope
SIF <- data$SIF
SSRT <- data$SSRT
sNR <- data$sNR
df <- gather(data, run, acc, run1:run8, factor_key = TRUE)
df$sNR <- as.factor(df$sNR)
#----------------
## ANOVA
contrasts(df$run) <- contr.poly
res.aov <- aov(acc ~ run + Error(sNR), data = df)
print(summary(res.aov, split = list(run = list(Linear = 1))))
## PLOT
p1 <- plotAcc(df, c("Runs"), "")
## OUTLIERS
sif_slope_outliers <- get_outliers(SIF, "SIF", "%", Slope, "Slope", "", disp = FALSE)
ssrt_slope_outliers <- get_outliers(SSRT, "SSRT", "%", Slope, "Slope", "", disp = FALSE)
## CORRELATIONS
# SIF/Slope
p2 <- plotCorrelation(SIF, Slope,"Inhibiting memories (SIF, %)","Classification accuracy slope",
doCorrelation(SIF, Slope, sif_slope_outliers),
1.8, 14, sif_slope_outliers, FALSE)
#SSRT/Slope
p3 <- plotCorrelation(SSRT, Slope,"Inhibiting actions (SSRT, ms)","Classification accuracy slope",
doCorrelation(SSRT, Slope, ssrt_slope_outliers),
1.8, 14, ssrt_slope_outliers, FALSE)
# plots <- grid.arrange(p1, p2, p3, nrow=2)
## Put all plots in a grid and display
plots <- grid.arrange(p1, p2, p3, nrow=3)
#plots <- grid.arrange(
# arrangeGrob(p1,
# top = "Accuracy of the action-stopping classifier to discriminate thought suppression",
# layout_matrix = rbind(c(1,1))),
# arrangeGrob(p2, p3,
# top = '\nRelationship between the accuracy declain and behavioural measures of inhibition\n',
# nrow = 1))
#
return(plots)
}
# ----------------------------------------------------------------------
# GET THE DATA
# ----------------------------------------------------------------------
ccdata <- read.csv('https://raw.githubusercontent.com/dcdace/tmp/master/data/perRunCrossClassification.csv')
ccdata[,5:12] <- ccdata[,5:12]*100
head(ccdata)
ccdata_rDLPFC <- subset(ccdata, ROI == 'rDLPFC')
ccdata_rVLPFC <- subset(ccdata, ROI == 'rVLPFC')
ccdata_HC <- subset(ccdata, ROI == 'HC')
ccdata_M1 <- subset(ccdata, ROI == 'M1')
options(repr.plot.width = 5, repr.plot.height = 10, repr.plot.res = 200)
rDLPFC_plots <-perRunResults(ccdata_rDLPFC)
#save
ggsave(file="FigureS4.pdf", rDLPFC_plots, width=5, height=10, device = cairo_pdf)
#cairo_pdf deals with font embedding
options(repr.plot.width = 5, repr.plot.height = 10, repr.plot.res = 200)
rVLPFC_plots <-perRunResults(ccdata_rVLPFC)
#save
ggsave(file="FigureS5.pdf", rVLPFC_plots, width=5, height=10, device = cairo_pdf)
#cairo_pdf deals with font embedding
options(repr.plot.width = 5, repr.plot.height = 10, repr.plot.res = 200)
HC_plots <- perRunResults(ccdata_HC)
options(repr.plot.width = 5, repr.plot.height = 10, repr.plot.res = 200)
M1_plots <- perRunResults(ccdata_M1)
# ----------------------------------------------------------------------
# GET THE DATA
# ----------------------------------------------------------------------
df_dcm <-
read.csv(
'https://raw.githubusercontent.com/dcdace/tmp/master/data/DCM_results.csv',
fileEncoding = "UTF-8-BOM"
)
df_dcm$Name <- factor(df_dcm$Name, levels = df_dcm$Name)
models <- subset(df_dcm, Type == 'm')
models$Name <- models$Nr
fDirection <- subset(df_dcm, Type == 'Direction')
fPathways <- subset(df_dcm, Type == 'Pathways')
fInteractions <- subset(df_dcm, Type == 'Interactions')
fTargets <- subset(df_dcm, Type == 'Targets')
# ----------------------------------------------------------------------
# PLOTTING FUNCTION
# ----------------------------------------------------------------------
plotDCM <- function(dataset, r, plottitle, ytitle) {
ggplot(dataset, aes(x = Name, y = EP, label = round(EP, r))) +
ggtitle(plottitle) +
geom_bar(
stat = "identity",
color = "black",
fill = "darkgrey",
alpha = 0.3,
lwd = 1) +
geom_text(data = dataset[dataset$EP > 0.03,], size = 5, nudge_y = 0.065) +
scale_y_continuous(ytitle,
limits = c(0, 1.1),
breaks = seq(0, 1, 0.2)) +
theme_classic(base_size = 20) +
theme(axis.line = element_line(size = 1),
axis.ticks = element_line(size = 1),
axis.text.x = element_text(angle = 35, hjust = 0.99),
axis.title.x = element_blank())
}
bp1 <- plotDCM(fDirection, 4, "Direction", "Exceedance probability")
bp2 <- plotDCM(fPathways, 4, "Pathways", "") + theme(axis.title.y = element_blank())
bp3 <- plotDCM(fInteractions, 2, "Interactions", "") +
scale_x_discrete(labels =
c("None",
expression("DLPFC" %->% "VLPFC"),
expression("DLPFC" %<-% "VLPFC"),
expression("DLPFC" %<->% "VLPFC"))) +
theme(axis.title.y = element_blank())
bp4 <- plotDCM(fTargets, 2, "Targets", "") + theme(axis.title.y = element_blank())
bp5 <- plotDCM(models, 2, "All models", "") +
scale_x_continuous(breaks = seq(0, 72, 4)) +
theme(axis.title.y = element_blank(), axis.text.x = element_text(angle = 0, hjust = 0.5))
# ==================================================================
# PLOT
# ==================================================================
options(repr.plot.width = 20, repr.plot.height = 5, repr.plot.res = 200)
plot.dcm <-
plot_grid(bp1, bp2, bp3, bp4, bp5, nrow = 1, align = "h",
rel_widths = c(6/15,5.5/15,5.5/15,3.3/15,12/15))
# display
plot.dcm
#save
svg("Figure6b_f.svg")
plot.dcm
dev.off()
# reset to defult plot size
options(repr.plot.width = 7, repr.plot.height = 7)